#
# pcfcross3D.R
#
# Cross-type pair correlation function for 3D multitype point patterns
#
# $Revision: 1.0 $
#

"pcfcross3D" <-
function(X, i, j, ...,
         rmax=NULL, nrval=128,
         correction=c("translation", "isotropic"),
         delta=NULL, adjust=1, biascorrect=TRUE,
         ratio=FALSE)
{
  # Verify input is a 3D multitype pattern
  if(!inherits(X, "pp3"))
    stop("X must be a pp3 object")
  if(!is.multitype.pp3(X))
    stop("X must be a multitype pp3 object (marks should be a factor with multiple levels)")

  if(is.NAobject(X)) return(NAobject("fv"))
  if(missing(correction))
    correction <- NULL

  # Extract marks and validate i, j
  marx <- marks(X)
  if(missing(i))
    i <- levels(marx)[1]
  if(missing(j))
    j <- levels(marx)[2]

  if(!(i %in% levels(marx)))
    stop(paste("Mark level", i, "not found in pattern"))
  if(!(j %in% levels(marx)))
    stop(paste("Mark level", j, "not found in pattern"))

  I <- (marx == i)
  J <- (marx == j)

  # Check if there are enough points of each type
  nI <- sum(I)
  nJ <- sum(J)
  if(nI < 2)
    stop(paste("Not enough points of type", i, "(", nI, "points found)"))
  if(nJ < 2)
    stop(paste("Not enough points of type", j, "(", nJ, "points found)"))

  Iname <- paste("points with mark i =", i)
  Jname <- paste("points with mark j =", j)

  # Call the underlying multitype function
  result <- pcfmulti3D(X, I, J, ...,
                      rmax=rmax, nrval=nrval,
                      correction=correction,
                      delta=delta, adjust=adjust, biascorrect=biascorrect,
                      Iname=Iname, Jname=Jname,
                      ratio=ratio)

  # Rebadge the output
  iname <- make.parseable(paste(i))
  jname <- make.parseable(paste(j))
  result <-
    rebadge.fv(result,
               substitute(g[i,j](r),
                          list(i=iname,j=jname)),
               c("g", paste0("list", paren(paste(iname, jname, sep=",")))),
               new.yexp=substitute(g[list(i,j)](r),
                                   list(i=iname,j=jname)))
  return(result)
}

pcf3engineCross <- function(x1, y1, z1, x2, y2, z2, box=c(0,1,0,1,0,1),
                           rmax=1, nrval=100, correction="translation",
                           delta=rmax/10)
{
  code <- switch(correction, translation=0, isotropic=1)
  res <- .C(SE_Rcallpcf3cross,
            as.double(x1), as.double(y1), as.double(z1), as.integer(length(x1)),
            as.double(x2), as.double(y2), as.double(z2), as.integer(length(x2)),
            as.double(box[1]), as.double(box[2]),
            as.double(box[3]), as.double(box[4]),
            as.double(box[5]), as.double(box[6]),
            as.double(0), as.double(rmax), as.integer(nrval),
            f=double(nrval),
            num=double(nrval),
            denom=double(nrval),
            as.integer(code), as.double(delta))
  return(list(f=res$f, num=res$num, denom=res$denom))
}

"pcfmulti3D" <- function(X, I, J, ...,
                         rmax=NULL, nrval=128,
                         correction=c("translation", "isotropic"),
                         delta=NULL, adjust=1, biascorrect=TRUE,
                         Iname="points satisfying condition I",
                         Jname="points satisfying condition J",
                         ratio=FALSE)
{
  verifyclass(X, "pp3")
  if(is.NAobject(X)) return(NAobject("fv"))

  correction.given <- !missing(correction) && !is.null(correction)
  if(is.null(correction))
    correction <- c("translation", "isotropic")
  correction <- pickoption("correction", correction,
                           c(translation="translation",
                             trans="translation",
                             isotropic="isotropic",
                             iso="isotropic",
                             best="isotropic"),
                           multi=TRUE)

  # Validate subset indices
  I <- ppsubset(X, I, "I")
  J <- ppsubset(X, J, "J")
  if(is.null(I) || is.null(J))
    stop("I and J must be valid subset indices")

  nI <- sum(I)
  nJ <- sum(J)
  if(nI == 0) stop(paste("There are no", Iname))
  if(nJ == 0) stop(paste("There are no", Jname))

  # Extract coordinates for subsets
  coo <- coords(X)
  cooI <- list(x=coo$x[I], y=coo$y[I], z=coo$z[I])
  cooJ <- list(x=coo$x[J], y=coo$y[J], z=coo$z[J])

  # Domain information
  B <- X$domain
  if(is.null(rmax))
    rmax <- diameter(B)/2
  r <- seq(from=0, to=rmax, length.out=nrval)

  # Intensity calculations
  volB <- volume(B)
  lambdaI <- nI/volB
  lambdaJ <- nJ/volB
  nIJ <- sum(I & J)
  samplesize <- npairs <- nI * nJ - nIJ
  lambdaIJarea <- npairs/volB

  # Bandwidth selection
  if(is.null(delta)) {
    delta <- adjust * 0.26/lambdaJ^(1/3)
  }

  # Bias correction
  if(biascorrect) {
    rondel <- r/delta
    biasbit <- ifelseAX(rondel > 1, 1, (3/4)*(rondel + 2/3 - (1/3)*rondel^3))
  }

  # Initialize output fv object
  g <- data.frame(r=r, theo=rep.int(1, length(r)))
  desc <- c("distance argument r", "theoretical Poisson %s")
  g <- fv(g, "r", quote(g[3](r)),
          "theo", , c(0,rmax/2),
          c("r", "{%s[%s]^{pois}}(r)"),
          desc, fname=c("g", "3"))

  # Extract domain bounds
  flatbox <- unlist(B[1:3])

  # Compute for each correction method
  if(any(correction %in% "translation")) {
    u <- pcf3engineCross(cooI$x, cooI$y, cooI$z,
                        cooJ$x, cooJ$y, cooJ$z,
                        flatbox,
                        rmax=rmax, nrval=nrval,
                        correction="translation", delta=delta)
    gt <- u$f
    if(biascorrect)
      gt <- gt/biasbit
    g <- bind.fv(g, data.frame(trans=gt),
                 "{hat(%s)[%s]^{trans}}(r)",
                 "translation-corrected estimate of %s",
                 "trans")

    # Store numerator and denominator for ratio if requested
    if(ratio) {
      attr(g, "trans_num") <- u$num
      attr(g, "trans_denom") <- u$denom
    }
  }

  if(any(correction %in% "isotropic")) {
    u <- pcf3engineCross(cooI$x, cooI$y, cooI$z,
                        cooJ$x, cooJ$y, cooJ$z,
                        flatbox,
                        rmax=rmax, nrval=nrval,
                        correction="isotropic", delta=delta)
    gi <- u$f
    if(biascorrect)
      gi <- gi/biasbit
    g <- bind.fv(g, data.frame(iso=gi),
                 "{hat(%s)[%s]^{iso}}(r)",
                 "isotropic-corrected estimate of %s",
                 "iso")

    # Store numerator and denominator for ratio if requested
    if(ratio) {
      attr(g, "iso_num") <- u$num
      attr(g, "iso_denom") <- u$denom
    }
  }

  # Set up output
  formula(g) <- . ~ r
  unitname(g) <- unitname(X)
  attr(g, "delta") <- delta

  # Handle ratio parameter
  if(ratio) {
    # Store numerator and denominator information as separate attributes
    # to avoid creating recursive ratio objects that cause C stack overflow
    # Keep the original attribute names to maintain compatibility
    # The attributes are already stored as trans_num/trans_denom and iso_num/iso_denom
  }

  return(g)
}